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ABSTRACT 


Several path following algorithms based on the combination of three smooth 
penalty functions, the quadratic penalty for equality constraints and the quadratic 
loss and log barrier for inequality constraints, their modern counterparts, augmented 
Lagrangian or multiplier methods, sequential quadratic programming, and predictor- 
corrector continuation are described. In the first phase of this methodology, one 
minimizes the unconstrained or linearly constrained penalty function or augmented 
Lagrangian. A homotopy path generated from the functions is then followed to 
optimality using efficient predictor-corrector continuation methods. The continuation 
steps are asymptotic to those taken by sequential quadratic programming which can 
be used in the final steps. Numerical test results show the method to be efficient, 
robust, and a competitive alternative to sequential quadratic programming. 
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and third authors were supported in part by the Air Force Office of Scientific 
Research through Grant # AFOSR-88-0059. 
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1. Introduction 


Path following algorithms for the solution of constrained optimization have been 
revitalized in recent years, due in no small part to the claims and success of the 
Karmarkar algorithm in linear programming. The ones presented here can be viewed 
as a combination of various elements and techniques in nonlinear programming: three 
smooth penalty functions (quadratic penalty for equality constraints, quadratic loss 
and log barrier for inequality constraints), their more modern counterparts, 
augmented Lagrangian or multiplier methods, sequential quadratic programming 
(Newton's method), and predictor-corrector continuation methods for efficient path 
following. The objective in this work then is to describe this class of algorithms and 
to present numerical evidence of their efficiency, robustness, and potential. 

One view (ref. 1) of these algorithms starts with the three aforementioned 
smooth penalty functions. One first performs an unconstrained or linearly constrained 
optimization of the penalty function. The minimizer then satisfies a set of first order 
necessary conditions (the gradient of the penalty function is zero when all constraints 
are incorporated into the penalty function) from which one can define an equivalent 
system of parameterized nonlinear equations. This system represents a perturbation of 
the Karush-Kuhn-Tucker first order necessary conditions, and the solution is followed 
to optimality using efficient predictor-corrector continuation methods. The simplest 
predictor-corrector steps are asymptotic to ' those taken in sequential quadratic 
programming, and thus the local convergence rates are the same as those of 
sequential quadratic programming. When shifts and weights are added to these three 
penalty functions and are adaptively chosen or updated during the optimization 
phase, one has the class of multiplier or augmented Lagrangian methods. 
Theoretically, one can expect a shorter path through the use of augmented 
Lagrangians (ref. 2, Theorem 12.2.1), which suggests that the use of these updates in 
the weights, scales, and shifts may be used to generate good paths to optimality. 

A different perspective of these algorithms evolves from sequential quadratic 
programming (SQP) itself. These SQP methods perform exceptionally well in 
minimizing function evaluations, but may be slow since the combinatorial complexity 
of the inequality constraints is reflected in the subproblems at each step. 
Furthermore, they are currently restricted to small to medium size problems with 
promise for large scale applications (refs. 2 and 3). Augmented Lagrangian methods, 
on the other hand, are currently used effectively for large scale problems with 
structured sparsity arising, for example, from discretized partial differential equations. 
Also, augmented Lagrangians are often used as merit functions for globalizing 
sequential quadratic programming. However, the minimizer of the augmented 
Lagrangian at any given stage is not a solution of the original problem, and thus the 
homotopy between the minimizer of the augmented Lagrangian and sequential 
quadratic programming may be viewed as an intermediate globalization technique. 
We find this to be very efficient. 

In the sections to follow, we briefly outline the methodology and present in 
section seven the results of our preliminary numerical testing. 
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2. Background and General Results 

For expendiency and convenience of presentation, we present in this 
section the first order necesssary and second order sufficient conditions for 
the mathematical programming problem 

(2.1) Min { f(x) | h(x) = 0, g(x) > 0} 

where f: IR n -> h: IR n -* 1R^ and g: IR n -+ IR^ are assumed to be twice continuously 
differentiable in an open set fi containing the feasible region 91 = { x | h(x) = 0, 
g(x) > 0 }. The Karush-Kuhn-Tucker first order necessary conditions state that in 
the presence of a constraint qualification, there exist multiplier vectors A and p such 
that 


( 2 . 2 ) 


F(x,A,/i) = 


V£(x, \,p) 
h( x ) 


[Mg(x) 

M > 0, g(x) > 0, 


= 0 


where C = f - 


h T A - g T /*, 


M 


= diag(/i), n = (Mp.-.^J 1 , A 



THEOREM 2.1 (ref. 4). Let (x be a solution of F(x, \,p) = 0. Assume f g 

and h are twice continuously differentiable in a neighborhood of x Q and define two 

index sets A and A and a corresponding tangent space T by 
A = {i: 1 < i < p, gfij = 0), A = {* G A: p° + 0 }, 

T = {y e K n : DJifxJy = 0, Dgfrjy = 0 (i G A)}. 

Then a necessary and sufficient condition that D^ x ^ ^( Z 0 >^ 0 >1 1 0 ) & e nonsingular is 
that each of the following three conditions hold: 

(a) A = A; 

(b) {{V ^-(x U {V^/V)p is a linearly independent collection of 

q + vectors where \A\ denotes the cardinality of A; 

2 

(c) the Hessian of the Lagrangian V X C is nonsingular on the tangent space T 

at (VW- 

Furthermore, if (b) remains valid, p Q > 0, gfx^ > 0, and (a) and (c) are 
replaced by 

(a') p°i > 0 for all i G A, 
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2 

(d) the Hessian of the Lagrangian ^ Jl is positive definite on the tangent 

space T at (x & \ & pj, 

then D( x ^ p)^( z o’^o’^ () ?s nonsingular and x Q is a local minimum of the nonlinear 
programming problem (2.1). 

We note that condition (a) is a complementarity condition, (b) is the linear 
independence constraint qualification, and (c) is a second-order condition. 


3. Newton's Method for Nonlinear Programming 

The objective of this section is to briefly review Newton's method 
(Newton-Lagrange or sequential quadratic programming) for the general nonlinear 
programming problem (2.1). The methods which we later present are asymptotic to 
Newton's method and, in fact, may be used in the last few iterations. The starting 
point for Newton's method is the first order necessary condition as stated in equation 
( 2 . 2 ). 


The motivation for Newton's method is to form a linear model suggested by 
Taylor's formula about the current point (x,A,//) and then solve this moael for the 
linear increments (Ax,AA,A^). The difference between a nonlinear equations 
approach and that of constrained optimization is that one keeps a quadratic term in 
the complementarity condition: 


(3.1) F(x+Ax,A+AA,/H- Ap) = 


V 2 LAx - D h T Ax - D g T Ax + VL 

A A 

D hAx + h 

A 

(M + AM)(D x gAx + g) 


+ 0(A 2 ) 


o 

where 0(A ) refers to those increments (Ax,AA,A//) which are second order and 
higher. Note that the quadratic increment AMD x gAx has been retained, but the 

addition of such higher order terms will not affect the quadratic convergence of 
Newton's method. Setting this bracketed quantity to zero yields 

V 2 LAx - D x h T Ax - D x g T Ax = -VL 

D x hAx + h = 0 

(M + AM)(D x gAx + g) = 0, 

which, if p + A// > 0 is imposed, represent the first order optimality conditions for 
the quadratic programming subproblem 

MIN f + Vf T Ax + (1/2 )Ax T V 2 LAx 
( 3.2) ST D x hAx + h = 0 

D x gAx + g > 0. 
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Second order information for both the constraints and the objective function are built 

2 

into the quadratic programming subproblem through the appearence of V L. The 
generic Newton's method with no safeguards is as follows: 

ALGORITHM [NEWTON'S METHOD FOR NONLINEAR PROGRAMMING] 

Initialize x = x . A = A . and u = « 
o’ o’ r r o 

For k = 0,1,..., until satisfied, do (a), (b), (c) 

(a) Compute f = f(x k J, h = h(x k ), g = g(x k ), Vf = Vf(x fc ), 

D x h = D x h(x k ), D x g = D x g(x k ), V 2 L = V 2 L(x k ,A k ,^), 

(b) Solve (3.2) for the correction Ax in x k and the multipliers 
/ 1 k + Ap > 0 and A k + AA, 

(c) Update: x k+1 = x k + Ax, A k+1 = A k + AA, ^ k+1 = p^+Ap. 

The quadratic convergence of Newton's method for nonlinear equations is preserved 
for the nonlinear programming problem under conditions similar to those in nonlinear 
equations: 

Theorem 3.1 [Convergence of Newton’s Method]. Let x be a local solution of the 

nonlinear programming problem, assume that f g, and h are cf functions whose 

second derivatives are Lipschitz continuous in a neighborhood of x, and suppose that 
the linear independence constraint qualification (condition (b) in Theorem 2.1) is 

satisfied. Then there exist multipliers A and p >0 for which the Karush- Kuhn- Tucker 
conditions (2.2) are valid. If conditions (a') and (b') in Theorem 2.1 are satisfied, 

then there is a neighborhood N of (x,X,p) such that if (x Q ,X o ,pJ is in N, then the 
iterates (x^X^pjJ are defined, remain in N and converge quadratically to (x,X,p). 

As with Newton's method for nonlinear equations, there are the questions of 
linear algebra, updating techniques for quasi-Newton steps, and globalization methods. 
These questions are discussed, for example, in the review paper of Stoer (ref. 3), the 
book of Fletcher (ref. 2), and the references therein. Two commonly used merit 
functions for globalization are the Lj penalty function 

P = Vt + (1/r) S |hj(x)| + (1/r) £ gj (xf 
and the augmented Lagrangian 

L a (x,A,//,r) = f + (l/2r) E (h.(x) - rAj) 2 + (l/2r) £ [(gj(x) - r/i.) _ f 
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where gj(x) = Min { gj(x), 0 }, r is the penalty parameter, and A and fi are the 

shifts or approximate multipliers. The use of the penalty function as a merit 

function for sequential quadratic programming has proved to be highly successful for 
small to medium problems. Augmented Lagrangians on the other hand have as their 
domain of application large scale problems such as those that appear in discretized 
partial differential equations; however, the function evaluation count on small to 
medium problems is not nearly as favorable as that of sequential quadratic 
programming (ref. 5). On the other hand, sequential quadratic programming may be 
slow in comparison to augmented Lagrangian methods, primarily due to the 
combinatorial complexity of the inequality constraints which persists in the quadratic 
subprograms. In this work we combine both approaches by using a few steps of 
augmented Lagrangians followed by a homotopy phase, and then sequential quadratic 
programming. 

4. Homotopy Methods for Constrained Optimization 

The idea of a homotopy method is to embed a difficult problem into a 
parameterized set of problems such that at one parameter value the problem is 
"easy" to solve and at another value one recovers the "difficult" problem. One then 
continues the solution of this parameterized system from the easy problem to the 
desired one. For nonlinear equations F(x) = 0, two commonly used homotopies are 
G(x,t) = tF(x) + (l-t)(F(x) - F(a)) (the global Newton homotopy) and 

G(t,x) = tF(x) + (l-tj(x-a) (the regularizing homotopy). Indeed, at t = 0 these 
are easy to solve, while at t = 1 one recovers the original problem. These 
homotopy methods tend to be quite robust, but currently are not as efficient as the 
use of merit functions with a modified or quasi-Newton method. The homotopy 
methods discussed here are generally very efficient, but the the easy problem is not as 
easy as the above ones for nonlinear equations in that the "easy" problem requires 
the solution of an unconstrained or linearly constrained optimization problem. 

To illustrate the idea and for later numerical comparisons, we first consider the 
mixed quadratic penalty-log barrier function 

(4.1) P(x,r) = f + (l/2r) E h?(x) - r E ln(gj(x)) 
or in the more general form 

(4.2) P(x,i/,A,^,<7,a>,£,r) = iA + (l/2r) E a.(7.h.(x) - rtr^Aj) 2 

- r E a/j ln(gj(x) + rtfj) 


wherein weights <7j and oij, scales v and 7 j, and shifts A and have been 

introduced. The homotopy generated from this penalty function depends on the 
system parameters, which can be adaptively chosen during the optimization phase. 


To explain how one can derive a homotopy, we consider the simpler form 
(4.1). At r = r , a minimizer, x , of this penalty function satisfies 
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VP = Vf + E Vhj(h,/r) - E (r/ gj )V gj = 0, 

which along with the definitions Aj = — h-/ r and ^ = r/gj yield an equivalent 
system of parameterized nonlinear equations 

VL = 0, L = f - h T A - g T n 

(4.3) h + rA = 0, 

Mg - re = 0, M = diag(p), e = (1,...,1) T . 

A solution to this parameterized system at r = r Q is given by A = -h(x Q )/r and 

fi = r/g(x ), componentwise. Furthermore, these equations represent the first order 

necessary conditions at r = 0 since /i(r) > 0 for r > 0. Once the optimization 
phase is complete, continuation techniques can be used to track the solution to 
optimality at r = 0. Further discussion of this homotopy can be found in the work 
of Poore and Al-Hassan (ref. 1). 

Another homotopy can be based on the quadratic penalty-loss function 

(4.4) P(x,r) = f + (l/2r) 2 h?(x) + (l/2r) 2 (g](x)] 2 
or more generally, the augmented Lagrangian function, 

(4.5) L a (x,!/,A,/i,(7,<M,Ar) = + (l/2t) 2 ^(^.(x) - r^'A,) 2 

+ (l/2r) E ^[(jSjgjfx) - r tfjif ) 2 

where gT(x) = Min { g.(x), 0 }, weights u, a- and <j)-, shifts Aj and jj... and scales 7 j 
and 0- have been introduced. These parameters can again be adaptively chosen 
during the optimization phase. At a minimizer of the augmented Lagrangian L & 

(4.6) VL a = M + (1/r) E Vf'jjh.K.r.-y.h^x) - rA,) 

+ (1/r) E V(7).g.)|(«»./).g.(x) - r/ij) - ] = 0. 

The definitions A = -(Srh - rA)/r and n = -((J>Bg - r jj) /r along with this equation 
yield the equivalent system 
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(4.7) 


VL = 0, 

rh + S -1 r(A - A) = 0 

/fc« + j - Aj) = o 


L = lA - (rh) T A - (Bg) 1 ^ 
S = diag(cr) 
if Ags(x) < fSfr. 

if ^jgj(x) > ^ 


where S = diag (a) and (j) = diag (<f>) represent the weights and T = diag (7) and 

B = diag (/?) the scales. The usual updates for the multipliers A and fi (refs. 2 

and 3) can be used in the optimization phase. Note that the use of different scales, 

weights, and updated multiplier approximations all change the homotopy path, and 
thus may be used to generate "good" paths. Furthermore, the ill-conditioning 

present in the penalty method is no longer present in the homotopies generated from 

these penalty functions. A final important modification to these homotopies is the 
normalization of the mulitpliers (ref. 1) to prevent multipliers tending to infinity, 
which happens generically when the linear independence constraint qualification is 
violated. 


5. Continuation Methods 

The system of parameterized equations posed in the previous section can be 
written as G(z,r)= 0 where r is the homotopy parameter and is arranged so that it goes 
from r Q > 0 to 0. The primary objective of this section is to briefly describe the 

methodology of traversing the path from r Q = r Q to optimality at r = 0. The idea 

is to generate a sequence of points {(zj, r j)}^_ 0 with r Q = r Q and r fl = 0. To get 

from (Zj,rj) to (z^pr^), one first predicts a new point (zP^pr?^) near the curve 

and then corrects back to the curve to obtain the desired ( z j_)_p r j + i)- Prediction is 

based on extrapolation of current and previous information about the solution. The 
extrapolation via polynomial interpolation of the solution values has been used for 
some time, but extrapolation of the tangents to the curve as is used here appears to 
be numerically more robust and efficient. A brief explanation of this methodology is 
given in the remainder of this section, but a more comprehensive explanation can be 
found in the works of Keller (ref. 6) and Shampine and Gordon (ref. 7). 

A formal differentiation of G(z,r) = 0 with respect to a third variable s yields 
Davidenko differential equation 


dz 
z- 3s 


D_G ^ + D 


dr 
r 3s 


0 , 


where s can be chosen to be arclength by adding the normalization 


dz m2 

3s 11 


+ 


dr 1 2 
3s I 


-1 = 0 . 
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dr 

Once an orientation, ie ^ is positive or negative, is known and if D z G is 
nonsingular, then one can write this system as the differential equation 


& - «w- 

Thus given a point (w^s^) on the curve, the diffential equation can be integrated to 
obtain 


w(s^ + As) = 



+ As 
f ( w) ds 


where f(w) denotes the tangent to the curve at the point w(s). If a polynomial 
Pk m of degree at most m is used to interpolate f at (wj c _j,Sj c _j) for j = o,...,m, the 

predicted solution is taken to be 

i Si + As 

w k+l = w k + As D{As), D(A S )=^J S * P km (s)d S . 


Given an error tolerance, once can vary the order of this formula to achieve the 
largest stepsize, As, possible. This method varies from the standard 
Adams-Bashforth technique in that the stepsize and order are varied at each step. 


Once a predicted point is obtained the correction back to the curve can be 
obtained in several wavs. Two popular ones are the vertical correction, wherein the 
system of equations G(z,r) = 0 is solved with r fixed at the predicted value, and the 
correction in a hyperplane orthogonal to the predictor direction. In this latter 
method one solves the augmented system of equations 

G(w) = 0 
N(w) = 0 


T P 

where N(w) = D (w - w^^) represents the plane orthogonal to the predictor 
direction and passing the predicted point w£ +1 and w = (z,r). 


6. Relation to Sequential Quadratic Programming 

The predictor-corrector method of the previous section gives steps toward the 
optimal solution and is, in fact, asymptotic to those obtained by sequential 

quadratic programming. The continuation phase may thus be viewed as a method 
for globalizing Newton's method. To explain the connection between the continuation 
phase with these two homotopies and sequential quadratic programming, we confine 
our attention in this section to the the homotopy generated by the quadratic 
penalty-log barrier function, ie. 
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VL = 0, 
h + rA = 0 
Mg - re = 0 


L = jA - h T A - g T n 


M = diag(/x). 

The result is that a Newton step as defined from the solution of (3.2) is asymptotic 

to a vertical correction plus an Euler prediction as r tends to 0. More precisely, if 

A = (Ax, AA, A n) denotes the Newton step as defined by the quadratic 

programming subproblems (3.2), Aj = (Axp AAp A/^), the Euler predictor or 

tangent to the curve with Ar = -r, and Ag = (Ax2 ? AA2, A^) the vertical 

correction with Ar = 0, then 

A = A^ + A2 + 0(r|| AA2II, ||A||^) as r -* 0 and A -* 0. 

A similar result applies for the quadratic penalty-loss function with the slight 
modification 


A — Aj + A2 + 0 (r|| A A25 A/^2 II ) as r -> 0. 

This suggests that as soon as the predicted value reaches r = 0, one could just as 
effectively switch to sequential quadratic programming without a globalizer. 

7. Numerical Examples 

In this section we consider two approaches to constrained optimization based on 
the quadratic penalty-log barrier function and the augmented Lagrangian or shifted 
quadratic penalty-quadratic loss function. For the former we have previously 
compared the numerics and briefly summarize some of their properties (ref. 1). For 
the augmented Lagrangian approach we present some of our recent testing on some 
nontrivial test problems (ref. 5) to demonstrate the robustness, efficiency, and 
potential for the methodology. 

For the quadratic penalty-log barrier function we first use the loss function 

(g ) g to generate a point x at which g(x) > 0 or is at least close to feasible 

region { x : g(x) > 0 } and then define a 6 so that g(x) + S > 0. Then we use a 
quasi-Newton with a BFGS update to minimize the penalty function P 

function P(x,r) = f(x) + h^(x)h(x)/(2r) - rS ln(gj(x) + r^/r°) at some value of 

the penalty parameter, say r , at which the problem is reasonably well conditioned. 

A quadratic-cubic line search and a Armijo stopping criterion (ref. 8), modified to 

maintain feasibility (g(x) 4- S > 0) has been used to globalize the quasi-Newton 
method. Once the minimization problem is solved, continuation techniques are used 
to track the solution to optimality at r = 0. The initial value of r Q = .1 has been 

used in the numerical experiments reported in the table below under the heading 
PENCON, but scaling has not been used. Additional information can be found in 
(ref. 1). 
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For the quadratic penalty-loss function, scaling and adaptive choices of the 
weights and scales have been used. For this penalty function, one does not need an 
initial feasible point for the inequality constraints. Again the BFGS update has been 
used, but the line search has been modeled after that of Fletcher (ref. 2). 

To get some estimation of the relative performance of this algorithm, we have 
solved several test problems from the book by W. Hock and K. Schittkowski (ref. 5) 
and give a comparative summary of the number of function evaluations in the table 
below. The function evaluation for the codes other than PENCON and LOSSCON 
are taken from (ref. 5). Consistent with those function evaluation counts, we count 
the evaluation of a p dimensional vector as p function evaluations; however, we do 
not count upper and lower bounds on variable since they are handled directly in the 
code and gradient evaluations of linear constraints are counted only once. The 
approximation of the Hessian of the Lagrangian in the continuation phase is based on 
finite differences (ref. 8). 


CODE 

AUTHOR 

METHOD 

VF02AD 

Powell 

Quadratic Approximation 

OPRQP 

Bartholomew-Biggs 

Quadratic Approximation 

GRGA 

Abadie 

Generalized Reduced Gradient 

VFOlA 

Fletcher 

Multiplier 

FUNMIN 

Kraft 

Multiplier 

FMIN 

Kraft, Lootsma 

Penalty 

PENCON 

Al-Hassan, Poore 

Penalty-Continuation 

LOSSCON 

Lundberg, Poore, Yang 

Loss Function-Continuation 


CODE: VF02AD OPRQP GRGA VFOlA FUNMIN FMIN PENCUN LUSSCON 

PROB. NO. 


5 

16 

16 

86 

32 

38 

200 

26 

16 

10 

48 

126 

678 

280 

554 

687 

88 

63 

12 

48 

132 

277 

300 

492 

306 

130 

68 

13 

180 

300 

192 

565 

928 

4,178 

269 

209 

14 

36 

126 

108 

192 

726 

838 

107 

91 

15 

30 

165 

508 

377 

496 

2,464 


113 

19 

78 

1,785 

314 

838 

3,339 

661 

423 

304 

20 

160 

200 

102 

383 

728 

4,094 

282 

104 

29 

52 

206 

646 

421 

482 

414 

155 

162 

118 

** 

1,800 c 

J , 857 

** 

** 

** 

2,804 



FUNCTION EVALUATION COUNT ^indicates failure 


Comparisons from this table illustrate that these penalty-continuation algorithms 
currently come in second and occasionally first and third; however, the methodology 
is quite robust, primarily because of the robustness of penalty paths leading to 
optimality and the robustness of the continuation methodology. The answers given 
by PENCON and LOSSCON have approximately twelve digits of accuracy on a 14+ 
digit CDC machine, whereas the remaining answers listed above are computed on a 
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10 digit machine (ref. 5) and those for VF02AD are generally to low accuracy. 
Furthermore, the quadratic penalty-loss function tends to perform a little better than 
the quadratic penalty - log barrier function with the current implementations. 
Scaling yields a significant improvement in some problems by reducing the number of 
function evaluations required in the unconstrained optimization phase and can reduce 
the number of steps taken in the continuation phase. The use of multiplier updates 
can also significantly reduce the length of the homotopy path but currently requires 
more function evaluations due to the reoptimization required after an update. These 
updates, however, may be useful for the large scale problems. 

8. Summary and Conclusion 

One of the objectives in this work has been to examine homotopy methods 
generated by three smooth penalty functions, the quadratic penalty for equality 
constraints and the log barrier and quadratic loss function for inequality constraints, 
and their modern counterparts, augmented Lagrangian or multiplier methods. We 
have shown that these methods are asymptotic to sequential quadratic programming 
and are competitive with these methods. However, considerable work will be 

required on large scale problems to assess their full potential. 

The robustness and efficiency of this class of algorithms based on these smooth 
penalty functions, augmented Lagrangians, the derived homotopy, and 

predictor-corrector continuation techniques have been illustrated in Section 7. The 
methodology shows considerable promise and potential for solving constrained 
optimization problems; however, as with any method a word of caution is 
appropriate. One can construct simple examples illustrating the following situations 
for penalty paths: given e > 0 a penalty path may exist only for the penalty 

parameter r > e or only for r < e, may not exist at all, may diverge, or may exist 

for r > 0 but the limit point at r = 0 is not a local minimum of the original 
problem. (When this last situation occurs, the Karush-Kuhn-Tucker equations have 
a singularity.) In spite of these examples penalty path following is a mathematically 
robust method of solving constrained optimization problems as is illustrated by the 
examples in Hock and Schittkowski (ref. 5), and these homotopy methods appear to 
make them efficient and a competitive alternative to the early stages of sequential 
quadratic programming. 
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